#ifndef PCPS_RANDOM_HEADER
#define PCPS_RANDOM_HEADER

#include <src/pcps.h>
#include <src/arena.h>

namespace pcps {

  //--------------------------------------------------------------------------------
  // pcps::lcg_int_t -- The integer type used in random number generators.
  //--------------------------------------------------------------------------------
  
  typedef unsigned int lcg_int_t;

  //--------------------------------------------------------------------------------
  // pcps::LinearCongruentialGenerator -- Class for a linear congruential
  //                                      random number generator.
  //
  //                                      Next number x = ( a * x + c ) mod m
  //--------------------------------------------------------------------------------

  class LinearCongruentialGenerator {

    private:

      // private data
      const lcg_int_t _m;
      const lcg_int_t _a;
      const lcg_int_t _c;
      lcg_int_t _x;

      // disable default constructor, copy constructor, and assignment operator
      LinearCongruentialGenerator();
      LinearCongruentialGenerator(const LinearCongruentialGenerator &);
      LinearCongruentialGenerator &operator=(const LinearCongruentialGenerator &);

    public:

      // constructor
      LinearCongruentialGenerator(const lcg_int_t m, const lcg_int_t a, const lcg_int_t c, const lcg_int_t seed)
        :  _m(m),
           _a(a),
           _c(c),
           _x(seed)
      {
        if (sizeof(lcg_int_t) != 4)
          throw pcps::Exception("pcps::LinearCongruentialGenerator requires lcg_int_t to be 32-bit");
      }

      // function to return the maximum possible value
      lcg_int_t max() const { return _m - 1; }

      // function to return the minimum possible value
      lcg_int_t min() const { return 0; }

      // function to return a random number
      lcg_int_t operator()() {
        const lcg_int_t retval = ( _a * _x + _c ) % _m;
        if (retval == _x)
          throw pcps::Exception("pcps::LinearCongruentialGenerator got stuck.  Try a different seed.");
        _x = retval;
        return retval;
      }

      // function to set the seed
      void set_seed(const lcg_int_t seed) { _x = seed; }

  };

  //--------------------------------------------------------------------------------
  // pcps::global_lcg -- A global linear congruential random number generator.
  //
  //                     Note that this object is NOT thread safe and should
  //                     be used only when the current thread has locked
  //                     pcps::rand_mutex().
  //--------------------------------------------------------------------------------

  extern pcps::LinearCongruentialGenerator global_lcg;

  //--------------------------------------------------------------------------------
  // pcps::LaggedFibonacci -- Class for a lagged fibonacci random number generator.
  //--------------------------------------------------------------------------------

  template <lcg_int_t p, lcg_int_t q> class LaggedFibonacci {

    private:

      // private data
      const lcg_int_t _word_length;
      const lcg_int_t _data_length;
      lcg_int_t _current_element;
      lcg_int_t * _data;

      // disable default constructor, copy constructor, and assignment operator
      LaggedFibonacci();
      LaggedFibonacci(const LaggedFibonacci &);
      LaggedFibonacci &operator=(const LaggedFibonacci &);

    public:

      // constructor
      LaggedFibonacci(pcps::Arena & arena)
        : _word_length(pcps::pow<lcg_int_t>(2,30)),
          _data_length(p+1),
          _current_element(p)
      {

        // check that values are sane
        if (sizeof(lcg_int_t) != 4)
          throw pcps::Exception("pcps::mc::LaggedFibonacci requires lcg_int_t to be 32-bit");
        if (q >= p)
          throw pcps::Exception("pcps::mc::LaggedFibonacci requires p > q");
        if (p <= 0)
          throw pcps::Exception("pcps::mc::LaggedFibonacci requires p > 0");
        if (q <= 0)
          throw pcps::Exception("pcps::mc::LaggedFibonacci requires q > 0");

        // get space for data array from the allocator
        _data = arena.allocate_array<lcg_int_t>(p+1, 0);

        // initialize the data vector
        pthread_mutex_lock(pcps::rand_mutex());
        bool have_odd = false;
        for (lcg_int_t i = 0; i < p; i++) {
          _data[i] = (lcg_int_t)( 0.1 + std::fabs( double(pcps::global_lcg()) * double(_word_length) / double(pcps::global_lcg.max()) ) ) % _word_length;
          if (_data[i] % 2 == 1)
            have_odd = true;
        }
        if (!have_odd)
          _data[pcps::global_lcg() % p] += 1;
        pthread_mutex_unlock(pcps::rand_mutex());

      }

      // function to return the maximum possible value
      lcg_int_t max() const { return _word_length - 1; }

      // function to return the minimum possible value
      lcg_int_t min() const { return 0; }

      // function to return a random number
      lcg_int_t operator()() {
        const lcg_int_t retval = (   _data[ ( (_current_element + _data_length) - p) % _data_length ]
                                   + _data[ ( (_current_element + _data_length) - q) % _data_length ] ) % _word_length;
        _data[_current_element] = retval;
        _current_element++;
        if (_current_element == _data_length)
          _current_element = 0;
        return retval;
      }

  };

  //-------------------------------------------------------------------------------
  // pcps::random_number -- returns a random number in the appropriate type
  //-------------------------------------------------------------------------------

  // no specialization
  template <class SCALAR> inline SCALAR random_number() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::random_number() for type \"%s\"")
                             % typeid(SCALAR).name() ).str() );
  }
  // specialization for double
  template <> inline double random_number<double>() {
    pthread_mutex_lock(pcps::rand_mutex());
    double retval = double(pcps::global_lcg()) / double(pcps::global_lcg.max());
    pthread_mutex_unlock(pcps::rand_mutex());
    return retval;
  }
  // specialization for complex
  template <> inline std::complex<double> random_number< std::complex<double> >() {
    return std::complex<double>(pcps::random_number<double>(), pcps::random_number<double>());
  }

  //-------------------------------------------------------------------------------
  // pcps::uni_01 -- returns a random number in the range 0 to 1
  //-------------------------------------------------------------------------------

  template <class GENERATOR> inline double uni_01(GENERATOR & gen) {
    const double retval = std::fabs( double(gen()) / double(gen.max()) );
    assert( retval >= 0.0 && retval <= 1.0 );
    return retval;
  }

  //-------------------------------------------------------------------------------
  // pcps::set_seed -- sets the random seed on each process
  //-------------------------------------------------------------------------------

  void set_seed(int seed);

} // end namespace pcps

#endif
